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Abstract 

Bose-Einstein condensates made of ultracold trapped bosonic atoms have become a central venue 
in which interacting many-body quantum systems are studied. The ground state of a trapped Bose- 
Einstein condensate has been proven to be 100% condensed in the limit of infinite particle number 
and constant interaction parameter [Lieb and Seiringer, Phys. Rev. Lett. 88, 170409 (2002)]. 
The meaning of this result is that properties of the condensate, noticeably its energy and density, 
converge to those obtained by minimizing the Gross-Pitaevskii energy functional. This naturally 
raises the question whether correlations are of any importance in this limit. Here, we demonstrate 
both analytically and numerically that even in the infinite particle limit many-body correlations can 
lead to a substantial modification of the variance of any operator compared to that expected from 
the Gross-Pitaevskii result. The strong deviation of the variance stems from its explicit dependence 
on terms of the reduced two-body density matrix which otherwise do not contribute to the energy 
and density in this limit. This makes the variance a sensitive probe of many-body correlations 
even when the energy and density of the system have already converged to the Gross-Pitaevskii 
result. We use the center-of-mass position operator to exemplify this persistence of correlations. 
Implications of this many-body effect are discussed. 
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I. INTRODUCTION AND MOTIVATION 


All physical information on the ground state of a quantum system can be obtained from 
its wavefunction by applying various operators and calculating expectation values. The 
variance of an operator quantifies to what extent the system under investigation is in an 
eigenfunction or a superposition of eigenfunctions of the corresponding operator. In this 
sense it sets the quantum resolution by which the corresponding operator could be measured. 
The variance of single-particle wavefunctions is analyzed in standard quantum mechanics 
textbooks, see for example [1]. How the many-body character of the quantum ground state 
effects the variance and how this effect persists with increasing the particle number will be 
the main points addressed in this work. 

Over the past two decades, since they were first experimentally realized [2-4], Bose- 
Einstein condensates (BECs) made of ultracold trapped bosonic atoms have become a pop¬ 
ular ground to study interacting quantum systems, see the reviews [5, 6] and books [7-9]. 
There has been tremendous theoretical interest in BECs, and ample studies have been made 
to describe them using Gross-Pitaevskii, mean-field theory. 

Correlations in quantum systems are responsible for many interesting phenomena. Gen¬ 
erally, one quantifies correlation as the deviation from a reference mean-field picture. The 
broadly accepted paradigm is that, in the limit of large particle number, Gross-Pitaevskii 
theory properly describes BECs in the ground state. To be more precise, the ground state 
of trapped BECs has been rigorously proven by Lieb and Seiringer to be 100% condensed 
in the limit of infinite particle number and constant interaction parameter [10]. Specifically, 
the energy and density of the condensate were shown to converge to those obtained by min¬ 
imizing the Gross-Pitaevskii, mean-held energy functional. This makes the ground state 
of a BEC a suitable system in which one can study the effects of correlations on different 
quantum observables and whether such effects persist in the infinite particle number limit. 

In the present work we pose the question whether there are properties of BECs which, 
unlike the energy and the density, are not obtained in the limit of infinite particle number 
from the Gross-Pitaevskii, mean-held result. If there were such properties, it would imply 
that there are many-body correlations in BECs that are not washed out - and that can be 
accessed - in the limit of infinite particle number! We show both analytically and numerically 
that the answer to this question is positive and that the variance of any operator is such 
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a property. The variance can deviate strongly from the Gross-Pitaevskii outcome even in 
the limit of infinite particle number. This demonstrates that the variance of an operator 
constitutes a sensitive probe for many-body correlations which are not washed out in this 
limit. 

In Sec. II we demonstrate how the energy and the density converge to the Gross-Pitaevskii 
result while the variance retains the many-body information of the system. The results 
are presented and discussed in Sec. Ill, with a sneak preview in Fig. 1 below. We stress 
that the conclusion holds true in the infinite particle number limit which is an exciting, 
counterintuitive result having a number of immediate and significant consequences. Of 
course, these correlations exist for BECs with a finite number of particles which are routinely 
produced in the laboratory. Summary and outlook are put forward in Sec. IV. 


II. THEORETICAL FRAMEWORKS 


A. Trapped Bose-Einstein Condensates in the Infinite Particle Limit 


We begin with the many-body Hamiltonian of N interacting bosons in a trap, 

N 

H{ t u ...,t n ) = ^h(r i ) + ^A°IT( r i - r fc ). (1) 

j=l j<k 

Here h( r) is the one-body Hamiltonian, with kinetic energy and trap potential terms, and 
W(v\ — r 2 ) the boson-boson interaction with Ao its strength. 

Consider the ground state of the system of bosons, 

H( ri,..., r 7V )'P(r 1 ,..., r N ) = EV( iq,..., rjy), (2) 


where the wavefunction ^(iq,..., iqv) is normalized to one and E is the energy. We shall 
find it instructive to continue our discussion by employing reduced density matrices of 
T(ri, ...,r N ) [11-13]. 

The reduced one-body density matrix is given by 

P (1) 0r,r'i) 


N 


= I dr 2 ...dr A r^*(r / 1 ,r 2 ,...,r J v)^(ri,r 2 ,...,r iV ) = ( 3 ) 


3 

The quantities aq(r) are the so-called natural orbitals and n 3 their respective occupations 
which are used to define the degree of condensation in a system of interacting bosons [14]. 
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FIG. 1. (Color online) The variance of X of a weakly-interacting Bose-Einstein condensate held 
in a symmetric trap for different barrier heights. Results for N = 1000 (in green), 10 000 (in 
blue), 100 000 (in magenta), and 1000 000 (in red) bosons are shown. The interaction parameter 
is A = \q{N — 1) = 0.1. (a) Shown is the variance -^A^, (full curves) and the Gross-Pitaevskii 
variance A| GP in red (in black; dashed curve). Large differences arise from certain barrier height, 
indicating that the correlations term A| corre / a tj ons becomes dominant. All four curves for the dif¬ 
ferent numbers of particles lie atop each other, (b) The many-body energy and (c) the depletion are 
seen to approach and coincide with the Gross-Pitaevskii results, as is expected from the literature. 
Note the small values on the y -axes of panels (b) and (c)! In contrast, the variance converges to a 
value substantially different from the Gross-Pitaevskii results for not too shallow barriers. See the 
text for further discussion. The quantities shown are dimensionless. 
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It is useful, although by no means mandatory, to express in what follows quantities using 
the natural orbitals oq. The (diagonal of the) reduced two-body density matrix is given by 

P = I dr 3 ...dr N V*(r 1 ,r 2 ,...,r N )V(r 1 ,r 2 ,...,r N ) = 

= - 1) a i( r l)°p( r 2)«fc( r l)«</( r 2), (4) 

jpkq ' ' 

where pj v k q = {'&\bjW p bkb q \'&) and the creation and annihilation operators are associated with 
the natural orbitals. The energy of the ground state is expressed as follows: 

77 = / dri 77 • «K ri )M r i)«i( r i) + 

j L 

+ N(T- 1) ' ~^~2 —“ / dr 2 a *j( r i) a *p( r z)W(ri - r 2 )a fc (ri)a p (r 2 ) . (5) 

pkq ^ ^ 

Note that in Eqs. (3-5) one-body quantities are divided by N and two-body ones by N(N— 1). 
The purpose will be seen shortly. We can now review the properties of the ground state in 
the infinite particle number limit [10] with the help of these reduced density matrices of the 
system. 

What is known on properties of the ground state in the infinite particle number limit, 
i.e., when N —> oo and the interaction parameter A = Xq(N — 1) is kept constant? With 
respect to the reduced one-body density matrix we have [10]: 

lim P ^ ri) = 7> GP (r 1 )^p(r' 1 ). (6) 

JV-»oo iv 

Equation (6) means that, in the limit of infinite particle number, the system of bosons is 
100% condensed, i.e., linqv-Kjo ^ = 1. Side by side, limjv-Kx> «i( r ) = 4>Gp( r ), namely, the 
first natural orbital becomes the Gross-Pitaevskii orbital. The latter minimizes the energy 
functional, 

dri0cp( r i) + —— J dr 2 (j)* G p(r 2 )W(r 1 - r 2 )7> GP (r 2 ) <f> G p(ri), (7) 

where £qp is the Gross-Pitaevskii energy. See for further discussion below. The density of 
the system is the diagonal of the reduced one-body density matrix, p(r) = p (1 ^(r,r). Then, 

lim dd = \fc P (r)\ 2 (8) 

JV->oo iv 

follows straightforwardly from Eq. (6). 
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With respect to the reduced two-body density matrix we have [10]: 


lim 

TV—>• oo 


p {2) {Y l ,Y2 1 Vl,r 2 ) 

N(N — 1) 


|0 G p(r 1 )| 2 |0 G p(r 2 )| 2 . 


( 9 ) 


Equation (9) means that Iirrrv->oo jvpv-q = 

Finally, combining the above infinite particle number limit of the reduced one-body and 
two-body density matrices, recalling that A = X 0 (N — 1) is held fixed, the energy of the 
ground state (5) becomes 


lim ^ = e GP , 
iV—>• oo TV 


( 10 ) 


where £gpi the Gross-Pitaevskii energy, is given in Eq. (7). 

We observe that the analysis of the infinite particle number limit [10] makes use of 
the dominant term, p m i, in the expansion of the reduced two-body density matrix. We 
henceforth pose the question what, if any, would be the consequence of the other terms? 
We would like to stress that such terms do not at all impact the literature result [10] on 
the density, energy, and the 100% condensation of the ground state, which, in the limit 
of infinite particle number, are, of course, described by the Gross-Pitaevskii results. But, 
perhaps, there are other properties of BECs which are not obtained in the infinite particle 
number limit from the Gross-Pitaevskii result? 


B. The Many-Body Variance 

Consider the ground state of N interacting bosons in a trap described by the wavefunc- 
tion ^(iq,..., rjv)- Let A(r) be an hermitian operator. In order to make our presentation 
more concrete, we choose A(r) = x, namely the position operator in the x direction. Other 
operators representing different quantities, for instance, the momentum operator, are infor¬ 
mative as well and the derivations below can be easily adapted. The choice of the position 
operator will also allow us for some immediate conclusions. 

We begin with the center-of-mass position operator, 

N 

x = Y.iy (n) 

3 =1 

The operator X is referred to as a one-body operator acting in the space of N particles. 
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A straightforward calculation gives the average of the center-of-mass position operator in 
the ground state, -h(T|A'|T) = f dr^-x, which is seen to be directly related to the density 
of the system. Using Eq. (8) we can readily obtain that 

= / dr I^Gp(r)| 2 ^- (12) 


Equation (12) states that, in the infinite particle number limit, the outcome of an average 
of a one-body operator in the ground state is given by the Gross-Pitaevskii result. 

But what about the variance of a one-body operator? To compute the variance we need 
the square of the center-of-mass position operator which is given by 

N 

X 2 = ^ Xj + ^2 (13) 

j =1 j<k 

and seen to be comprised of one-body and two-body operators. The presence of the two- 
body operator in (13) might suggest that the variance can capture many-body correlations 
in the interacting system. Recall, however, that the Hamiltonian (1) also contains a two- 
body operator, the boson-boson interaction. Yet, no many-body correlations are captured 
in the limit of infinite particle number by the energy [see Eq. (10)], so what could be the 
difference here? 

Expressed in terms of the above quantities, the variance of the center-of-mass position 
operator of the system takes on the form 


-1a| = i (<*|x 2 |w) - ('flu®) 2 ) = 


= f *Ak 2 _ N 


N 


N(N — 1) 


(N - 1) / dr 2 a*(r 1 )a*Jr 2 ) x 1 x 2 a k (r 1 )a q (r 2 ). 


Equation (14) deserves a closer look. The variance consists of two kinds of terms, one- 
body terms and a two-body term. We will study the impact and interplay of the one-body 
and two-body terms. Crucially, the latter is enhancedby the factor (N — 1) which will lead 
to intriguing results in the limit of infinite particle number. We note that, in the absence of 
boson-boson interaction, the variance boils down to that of a single particle. 

We can now turn to the infinite particle number limit of Eq. (14) which gives the final 
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expression 


1' —A 2 — A 2 -i- A 2 

^x,GP ' Wr, correlations': 


a I,gp = J dr\(j) G p(r)\ 2 x 2 - j dr|c/> GP (r)| 2 x , 


A l,correlations = Y1 / dr 2 x x x 2 a fc (ri)a,(r 2 ). (15) 

In obtaining Eq. (15) we have made used of the relations lim^v-^oo«i(r) = 0Gp( r ) and 
linijv-Kx> ^(n-i) = 1 discussed above. We see that, in the infinite particle number limit, the 
variance can be expressed as the Gross-Pitaevskii variance, A 2 GP , plus a second term which 
we call the correlations term, ^correlations- This is already fundamentally different than the 
situations for the quantities described so far, namely, the energy, density, and condensation. 
The question that arises is what is the nature of the second term in Eq. (15)? A central 
matter to pay attention to is that the summation therein only includes terms where at least 
one of the indexes is a higher natural orbital (i.e., second or above). Hence, A| corre iations 
only contributes in the presence of occupation of natural orbitals other than the first natural 
orbital (the so-called condensed mode). To remind, the depletion of the BEG is given by the 
occupation of all but the first natural orbital divided by N. Clearly, in the limit of infinite 
particle number the depletion becomes infinitesimal. Still the question remains whether a 
tiny occupation of the higher natural orbitals can make a difference? Yes, it can. Does this 
contradict the literature results for the energy, density, and condensation [10]? No, it does 
not. Let us see how. 

The occupation of the first natural orbital is of order O(N), and the (sum of the) rest is 
of order 0(1). This is compatible with the literature result that 100% condensation and the 
Gross-Pitaevskii density are recovered in the limit of infinite particle number, see Eqs. (6) 
and (8). The first matrix element of the reduced two-body density matrix, p im , is of order 
0(N 2 ), and the (sum of the) rest is of order O(N). This is compatible with the literature 
result that the Gross-Pitaevskii energy is recovered in the infinite particle number limit, see 
Eq. (10). But, unlike the case of the energy which is scaled by the interaction parameter 
A = Ao (N — 1), the variance is enhanced by the factor (N — 1). This allows the variance 
to pick-up, in the limit of infinite particle number, information from the other terms of the 
reduced two-body density matrix that do not contribute to the above quantities! 

We begin with an analytical case. Consider interacting bosons in an harmonic trap, an 
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isotropic one without loss of generality. The one-body Hamiltonian is h(r) = — \ (j^ + + 

) + |^ 2 r 2 , where u is the frequency of the trap and h — m — 1. The details of the boson- 
boson interaction, namely its strength and shape, are not important for our discussion. It 
is well known that one can separate the center-of-mass coordinate, - 7 = r ji from the 
remaining, (N — 1) relative-motion coordinates. Consequently, the iV-boson wavefunction 
becomes a product of the center-of-mass wavefunction, (-) e 1 , and a relative- 

motion wavefunction of the respective (N — 1) coordinates. It goes without saying that, due 
to this separability, the variance takes on the simplest form, 



1 

2u’ 


VN, 


(16) 


i.e., irrespective to the number of particles N. Eq. (16) thus holds (also) in the infinite 
particle number limit. 

On the other hand, it is well known that a repulsive boson-boson interaction broadens the 
Gross-Pitaevskii orbital, 4>gp( r), in comparison with that of the non-interacting system [5-9]. 
Generally, the Gross-Pitaevskii orbital depends on the interaction parameter, A = X 0 (N— 1), 
and on its shape. Thus, we arrive at the observation that A| GP for repulsive bosons in an 
harmonic trap is always larger than unless the interaction is zero, and then the 

occupation of the higher natural orbitals is strictly zero. As a consequence of this analysis, 
we see that it is the second term in Eq. (15), which is non-vanishing and negative, that 
accounts for the always existing many-body correlations in the system leading to the value 
(16) of the variance. To remind, the variance (16) is independent of the interaction strength 
and equals that of a single particle in an harmonic potential [1]. Thus, the relation 


A 2 

x, cor relations 


2 u 


- A 2 

a i,GP> 


(17) 


which is easily computable for a general system of interacting bosons in an harmonic poten¬ 
tial, can be used to quantify these correlations in the limit of infinite particle number. 


III. NUMERICAL INVESTIGATIONS AND DISCUSSION 

The results of the previous Subsec. IIB, when put in the context of the infinite particle 
number limit, are somewhat surprising. We have seen in the analytical example why the 
variance is not given in the limit of infinite particle number by the Gross-Pitaevskii result. 
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FIG. 2. (Color online) The density of a Bose-Einstein condensate held in a symmetric trap. The 
barrier height is (a) A=3 and (b) A=12. The number of bosons is N = 100 000 and the interaction 
parameter is A = Xq(N — 1) = 0.1. The density (in red; full curve) coincide with the Gross- 
Pitaevskii limit (in black; dashed curve), as is expected from the literature, whereas the respective 
variances deviate substantially from each other, as seen in Fig. 1. See the text for more details. 
The quantities shown are dimensionless. 


What happens to the variance when the bosons are trapped in an unharmonic potential, 
where the center-of-mass is not separable? In this case an analytical treatment is generally 
excluded and one is opt to resort to numerical investigations. 

We would now like to investigate numerically the variance of the center-of-mass posi¬ 
tion operator X of a BEC held in an unharmonic trap potential. We need a suitable and 
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accurate nmny-body tool to arrive at detailed conclusions. Such a many-body tool is the 
multiconfigurational time-dependent Hartree for bosons (MCTDHB) method, that has been 
well documented [15-18], extensively applied [19-23], and benchmarked [24] in the literature. 
To obtain the many-body ground state we propagate the MCTDHB eqnations-of-motion in 
imaginary time [24-26]. For the computation details, please see the appendix. 

We begin with a symmetric system in one spatial dimension. Specifically, we consider 
N = 1000, 10 000, 100 000, and 1000 000 bosons in a double-well potential. The one-body 
Hamiltonian, taken from Ref. [13], is h(x) = 2 + ^ + Ae~ x2 . For A = 0 there is no 

barrier, the trap potential is harmonic, and the variance is \ for any number N of particles, 
see Eq. (16). We would like to study the effect of the barrier on the variance for up to 
one million trapped bosons, and see what can be concluded on the variance in the infinite 
particle number limit. The interaction is taken to be the standard contact interaction, 
W(x 1 — x 2 ) = S(xi — x 2 ), with the interaction parameter A = A 0 (A^ — 1) = 0.1 held fixed. 
This guarantees that for any number of particles N the Gross-Pitaevskii, mean-field solution 
is the same. The results are collected in Fig. 1. 

Fig. la shows the variance -Tand the Gross-Pitaevskii variance A| GP for N = 1 000, 
10 000, 100 000, and 1 000 000 bosons. We can see that for low barrier heights the two quan¬ 
tities essentially coincide, whereas substantial differences emerge for higher barrier heights. 
As the number of particles increases, the difference between the ground-state energy and 
the Gross-Pitaevskii energy decreases, see Fig. lb. The depletion of the BEG decreases 
similarly, see Fig. lc. Note the small values on the y -axes of Figs, lb and lc! Side by side, 
the density of the BEC and the Gross-Pitaevskii density coincide, see Fig. 2. The behav¬ 
ior of the energy, the depletion, and the density seen in Figs. 1 and 2 are in line with the 
literature result on the infinite particle number limit. We can see that with increasing of 
the barrier height the depletion increases, and so does the deviation of the variance 7(7 A 2 _ 
from the Gross-Pitaevskii variance A| GP . This manifests the role of even a tiny amount 
of population outside the condensed mode leading to a strong impact on the variance, see 
Eq. (15). 

Returning to Fig. la, in contrast to the behavior of the depletion, the energy, and the 
density, the curves of the variance for different N are seen to lie atop each other. The 
latter indicates that the correlations term A| corrdations dominates the variance in the infinite 
particle number limit. An interesting point to take from Figs, lb and lc is how the limit of 
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infinite particle number is approached with increasing N. The respective curves are seen to 
be shifted vertically from each other. This indicates that, whereas the depletion decreases 
in the infinite particle number limit, the total number of particles residing in the higher 
natural orbitals is constant. When enhanced by the number of particles (N — 1) in the 
reduced two-body density matrix, this leads to the macroscopic difference of the variance 
from the Gross-Pitaevskii result, see Fig. la. 

To shed more light on the physical meaning of the variance A A'h we examine the densities 
shown in Fig. 2 more closely. For barrier height A — 3 the trap is essentially a single well. 
Accordingly, there is a single density peek, see Fig. 2a. The depletion for N = 100 000 
bosons is just above 10 -9 and consequently, the correlations part A% correlati(ms is very small. 
In this case, the variance essentially equals the Gross-Pitaevskii result, and both match the 
full-width-half-maximum of the density (~ 2y/r A 2 ^). When the barrier height is A = 12, 
the density is separated into two parts, see Fig. 2b, yet the system is essentially condensed, 
with depletion of slightly above 10~ 5 . This amounts to about 1 atom outside the condensed 
mode (the first natural orbital)! But that is enough, due to the enhancement by (N — 1) 
in the correlations term A? correlati(ms , the variance now differs substantially from the Gross- 
Pitaevskii limit. Interestingly, the variance matches the the full-width-half-maximum of each 
density peek (~ 2y / ^AV), whereas the Gross-Pitaevskii limit corresponds to the distance 
between the density peeks (~ 2^j A^ GP ). This is quite an effect of a single excited atom. 
And it persists in the limit of infinite particle number, because there will always be that 
excited atom, enhanced by the (N — 1) factor! 

Finally, the results plotted in Figs. 1 and 2 are obtained at the level of M = 2 orbitals 
in the MCTDHB method and found to be converged. For details of the convergence of the 
variance with increasing the number of orbitals M please see the appendix. 

We have so far treated the variance of X in a symmetric trap; Whether the harmonic 
trap, which admits an analytical treatment, or a double-well trap, which requires a numeri¬ 
cal investigation. When there is spatial symmetry in the system the natural orbitals possess 
symmetry as well, e.g., they are even and odd functions in case of a one-dimensional sym¬ 
metric trap. Consequently, there are terms of the reduced two-body density matrix that 
vanish due to the spatial symmetry of the system. 

As we have seen above, the difference in the variance in the infinite particle number 
limit from the Gross-Pitaevskii result has to do with excitations of atoms to higher natural 
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FIG. 3. (Color online) Same as Fig. 1 but for an asymmetric trap. Results for N = 100000 and 
the two interaction parameters A = Xq(N — 1) = 0.5 and 1 are shown. Note the small amount 
of depletion. The variance behave qualitatively the same as in the symmetric trap, compare to 
Fig. 1, and deviates substantially from the Gross-Pitaevskii limit for not too shallow barriers. The 
quantities shown are dimensionless. 

orbitals. No argumentation based on spatial (or other) symmetry was used, hence we do not 
expect this many-body effect to depend on symmetry. Still, the interplay of the different 
terms in Eq. (15) varies in the asymmetric case which warrants an example. 

The one-body Hamiltonian now takes the from h(x ) = — + ^ + Ae ~^ x ~ om ^ 2 / 8 . The 

interaction is taken again to be W(x i — x 2 ) = 5(xi — x 2 ) with the interaction parameters 
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FIG. 4. (Color online) Same as Fig. 2 but for an asymmetric trap. The number of bosons is 
N = 100 000 and the interaction parameter is A = Xo{N — 1) = 0.5. The density and the Gross- 
Pitaevskii result coincide. Compare to the respective variances in Fig. 3. The quantities shown are 
dimensionless. 


A = 0.5 and 1. The number of bosons is N = 100 000. The results are collected in Figs. 3 
and 4 and show that lifting the spatial symmetry does not change the conclusion that the 
variance of X of a trapped BEC can differ substantially from the Gross-Pitaevskii limit. Note 
that, despite the apparent slight asymmetry of the potential, the density is quite imbalanced 
between the two wells, see Fig. 4b. 

As a final example, we consider a finite BEC made of N = 100 bosons in a three- 
dimensional, unharmonic trap. We would like to touch upon how unharmonicity and inter- 
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action conspire to vary the variance. 

The one-body Hamiltonian now takes on the from h( r) = —| + + + ^ (x 6 +y 6 + 

z 6 ). The potential trap is isotropic and cubic-shaped. The short-range repulsive interaction 

/s _(r^ —r 2 )^ 

between the bosons is modeled by a Gaussian function W{ri — r 2 ) = ^^ 2 p /2 e 572 , where 

cr = 0.25 and the interaction parameters are A = Ao(A — 1) = 0.1, 1, 10, and 100. The 
results are presented in Fig. 5. As the interaction grows, the many-body variance decreases 
whereas the Gross-Pitaevskii result increases, see Fig. 5a. Here we plot and compare the 
radial variances, namely, -A(A-l + A 2 _ + A|) and A| GP + A ^ GP + A 2 GP . Side by side, 
the depletion increases, see Fig. 5b. Yet, even for the largest interaction we consider, the 
depletion is 1CU 2 , i.e., about 1 atom is excited outside the first natural orbital. Still, there 
is nearly a factor of one half when comparing the variance with the Gross-Pitaevskii result. 
Interestingly, the density, computed at the level of M = 4 orbitals with MCTDHB, see 
Fig. 5c, and the Gross-Pitaevskii result, see Fig. 5d, hardly differ. These densities clearly 
illustrate that, using the variance as a sensitive probe of many-body correlations, there is 
much into BECs than meets the eye. This concludes our investigations. 


IV. SUMMARY AND OUTLOOK 

Many-body correlations are an important part of any trapped interacting Bose-Einstein 
condensate. Even in the limit of infinite particle number, where the energy and density 
have converged to their Gross-Pitaevskii counterparts, the effect of the correlations persists 
and can be directly observed in the variance of any operator. In this sense the variance 
of an operator serves as a highly sensitive probe of the many-body correlations in a BEC 
comprised of any particle number. 

Taking a BEC center-of-mass position as an example, we demonstrated in a variety of 
systems that the variance of an operator can substantially differ from that predicted by 
the Gross-Pitaevskii theory even when only one in a million particles are excited out of the 
condensate. All in all, the results and conclusions of this work are expected to open new 
perspectives for research of trapped BECs. 

Further investigations on the dependence of the variance on the interaction strength 
and shape, as well as its dependence on the topology of the trap, will shed new light on 
many-body correlations in BECs. Studies of the variance of different observables, such 
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3.2 



FIG. 5. (Color online) (a) The radial variance (in red; full curve) of a three-dimensional Bose- 
Einstein condensate made of N = 100 atoms in an cubic-shaped unharmonic trap as a function 
of the interaction parameter A = \q(N — 1). The Gross-Pitaevskii variance is shown in black 
(dashed curve), (b) The difference between the many-body and Gross-Pitaevskii energy (in blue; 
full curve) and the depletion (in magenta; dashed curve), (c) The many-body density and (d) the 
Gross-Pitaevskii density for A = 100. Note the depletion is less than 1%. Shown are contours for 
4^ = 0.02 and 0.001. The many-body density coincide with the Gross-Pitaevskii density, whereas 
the respective variances deviate substantially from each other. See text for further details. The 
quantities shown are dimensionless. 
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as the momentum and angular-momentum operators, will provide further information on 
many-body correlations of trapped BECs. Last but not least, it is inviting to extend the 
investigations of the present work to the dynamics of BECs, which could only offer more 
insight into the importance of many-body correlations [27]. 
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Appendix A: Details of the Numerical Computations 

The MCTDHB method is well documented in the literature, ft is based on the real-space 
Hamiltonian H(ri ,..., r^) and utilizes M time-dependent orbitals which are determined 
according to the variational principle. As the method is variational, propagation in imagi¬ 
nary time convergences in the limit M —» oo orbitals to the ground state of the many-boson 
problem [24, 25]. On the other end, the Gross-Pitaevskii, mean-held theory is recovered for 
M =1. In turn, utilizing M time-dependent orbitals allows one to obtain accurate numerical 
results with substantially less numerical resources than a corresponding computation with 
M fixed orbitals would need. The accuracy of the method has recently been demonstrated 
[24], We use the implementation in [28]. 

In the present work the MCTDHB method is employed to compute the ground state in 
a one-dimensional double well (201 harmonic-oscillator discrete-variable-representation grid 
points in a box of size [-20,20]) and in a three-dimensional cubic-shaped unharmonic trap 
(64 3 grid points in a cube of size [-6,6] x [-6,6] x [-6,6] combined with a fast-Fourier-transform 
representation). As a concrete example and without loss of generality, convergence with 
increasing M of the variance of X of a Bose-Einstein condensate held in the one-dimensional 
symmetric trap of Fig. 1 is demonstrated in Fig. 6. 
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Barrier Height Barrier Height 


FIG. 6. (Color online) The variance of X of a Bose-Einstein condensate held in the symmetric trap 
of Fig. 1 for different barrier heights. Results for IV = 10 and 100 bosons for increasing number 
M of time-dependent orbitals are shown. The curves sit atop each other, and the variance seen to 
converge. The quantities shown are dimensionless. 
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